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Abstract. Phase separation of a two-dimensional van der Waals fluid subject to a grav- 
itational force is studied by numerical simulations based on lattice Boltzmann methods 
(LBM) implemented with a finite difference scheme. A growth exponent a. = 1 is mea- 
sured in the direction of the external force. 
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1 Introduction 

Phase ordering in fluids is an important process that still needs to be completely un- 
derstood in many cases of practical relevance. When a fluid is quenched from an initial 
disordered state into a regime of two-phase coexistence below the spinodal line, domains 
of the two phases are formed and grow with time. The typical size R of domains follows 
the power law R~t a with the growth exponent oc being universal in the sense that it does 
not depend on the microscopic details of the fluid, assuming only a few values related 
to the physical mechanism operating during phase separation [1]. Hydrodynamics is 
in general relevant and the coupling with the velocity field can change the value of the 
growth exponent a from that of purely diffusive growth [2,3]. 

In this paper we consider the ordering of a liquid-vapor system subject to an external 
field mimicking the effects of gravity. The role of gravity on phase ordering has been 
more studied in binary systems. In critical quenches, after an initial diffusive growth 
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with exponent 1/3, there is a viscous growth characterized by a. = 1 followed by an iner- 
tial regime with a=2/ 3 [4]. Gravity becomes relevant when heavy domains resting on top 
of light ones become gravitationally unstable, thus accelerating the domain growth [5]. 
This occurs at late stages making inertial growth difficult to observe. A theoretical anal- 
ysis neglecting hydrodynamical contributions suggests an exponent a z = l for the size 
of domains in the vertical direction [6]. There are few studies of phase separation for 
liquid-vapor systems. In two-dimensional simulations the values a = 1 /2 for high viscos- 
ity fluids and a = 2/3 for low- viscosity fluids have been found [2, 7]. We are not aware 
of simulations made on a liquid-vapor system subject to gravity where the growth expo- 
nent is measured. 

We address this problem by applying the lattice Boltzmann method (LBM) to simulate 
a van der Waals fluid described by the Navier-Stokes and the continuity equation. LBM 
have been proved successful in studying fluids with mesoscopic structures (liquid-vapor 
interfaces in our case) on large time scales, as it is needed for phase separation [2,4,8-12]. 
In our approach, the thermodynamic description is based on a free-energy functional 
where interfaces are described at a coarse-grained level. The free-energy interface cost is 
expressed, as usual in van der Waals-Landau models, in terms of gradients of the den- 
sity field. Locally, the fluid satisfies the van der Waals state equation. A finite difference 
version of LBM is implemented where the relationship c = 5s/5t among the lattice speed 
c and the space and time steps 5s and 5t does no longer hold, as in standard collision 
- streaming LBM [8-12]. The rejection of this condition has two advantages. First, this 
allows one to further consider multicomponent fluid systems where the masses of the 
component particles, as well as the lattice speeds, may be no longer identical [13, 14]. 
Second, higher order numerical schemes (including flux limiter schemes) may be consid- 
ered in order to reduce unphysical effects like the spurious velocity and the numerical 
viscosity [7, 13-18]. The use of high order numerical schemes in finite difference LBM 
helps further to improve the numerical stability and accuracy [7, 17] while providing a 
convenient alternative to interpolation supplemented LBM [19,20]. 

Our main results is that the sedimentation process induced by gravity is characterized 
by an exponent a. = 1 independently on the values of viscosity and gravity. 

The paper is organized as follows. Our LBM approach is described in Section 2; nu- 
merical results are shown in Section 3 and conclusions will be drawn in Section 4. 



2 Description of the model 

In this paper, we use the D2Q9 isothermal finite difference lattice Boltzmann model in 
two dimensions, which is well known in the literature [9-11,15,21]. This model relies on 
the following set of Af = 9 evolution equations for the non-dimensionalized distribution 
functions /i(r,f), i = 0,l,...,J\f— 1, defined in the nodes r= (x,y) of a lattice with A* x A y 
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nodes [7-12,18,21,22] 



fi(t,t+5t) = fiM-Stei-VfM-^lfiM-fffat)] 
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2 F-[ e! -u(r,f)]/; ei? (r,f). (2.1) 



In this non-dimensionalized model, the mass of particles equals 1. To reduce numerical 
errors, the second order Monitorized-Central-Difference flux limiter scheme [7,13,18,23, 
24] was used to calculate the space derivative V/)(r,f). 

The local values of the fluid quantities (particle number density n and velocity u) are 
derived from the distribution functions, as follows 

n(t,t) = M ^fi(T,t), (2.2) 

■(«>') = (2-3) 

The velocity vectors { e, } are given by 
e = 

e, = (cos^^,sin^±l)c for / = ! 4 (2.4) 

tt(2z-9) 

cos , sin cv2 for 1 = 5,. 



where c = c/cr = yO/x is a non-dimensionalized speed, 6=T /Tr is the non-dimensiona- 
lized system temperature, and # = 1/3. As discussed in [15,25], the following reference 
quantities «r = Na/V WC/ Tr = T c and cr = \Jk-gT c lm-R, where Na is Avogadro's number, 
V mc is the molar volume at the critical point of temperature T c , mg is the mass of the fluid 
particles and kg is Boltzmann's constant, may be used to get the non-dimensionalized 
values of the particle number density, temperature and speed, respectively. The system 
size is chosen as reference length Zr, the reference quantities £r and Ar for time and accel- 
eration follow from 

^ = 1 , ^ = 1. (2.5) 

Since we use finite difference schemes [7, 13, 15, 22, 23] to evolve the particle distribution 
functions according to Eq. (|2.1[) , the lattice spacing Ss and the time step Ss are no longer 
related to the lattice speed c as in standard lattice Boltzmann models [8-11], the temper- 
ature 6 is now a control parameter in our simulations. This feature of the finite difference 
approach allows us to change the system temperature 6 (and hence also the lattice speed 
c) while preserving the lattice spacing 5s and the system size in each direction (i.e., the 
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corresponding number of lattice nodes). In such models there is more freedom to choose 
the discrete velocity set, as done recently in a thermal model [26] where the possibility of 
having different sets of velocities allows to release the constraint of constant temperature. 

The equilibrium distribution functions f- q =f- q (r,t) that appear in the evolution equa- 
tion (|2.1|) of the D2Q9 model are expressed as series expansion up to second order with 
respect to the fluid velocity u [9-11,21] 



wm 



e,u (e,-u) 2 
1 + ^+ v 



where the weight coefficients are 
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(2.6) 



(2.7) 



The force F in Eq. (|2.1[) is introduced to recover the macroscopic equations of an 
isothermal van der Waals fluid subjected to the gravitational acceleration a. When us- 
ing the reference pressure p% = mn^c\, this force has the Cartesian components [7, 15,25, 
27-29] 

i 

(2.8) 



F« = -9«(p 
n 



' rfefl '-p Trafl ' s )+)ca a (V 2 n)+fl a 



where 



.ideal 



n6 



is the ideal gas pressure and 



.waals 



n6 



1 — bn 



-an 



(2.9) 



(2.10) 



is the van der Waals fluid pressure in non-dimensionalized form, k is a constant which 
controls the surface tension and a is the gravitational acceleration. The parameters a and 
b are given by 



7 U C 

8^ 

1 
3n~ c 



(2.11) 
(2.12) 



where n c is the particle density at the critical temperature 8 C = 1. In the following we will 
consider n c = 1 so the the van der Waals fluid pressure reads 



.waals 



3n6 _9 
3-n ~8' 



(2.13) 
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Note that the force (|2.8|) was already used in [7, 15, 18, 25] to investigate the dynamics 
and morphology of phase separation in liquid-vapour systems in the absence of gravity 
(a = 0). 

A Chapman-Enskog expansion shows that the continuity and Navier-Stokes equa- 
tions are recovered in the continuum limit: 

d t n+dp(nup) = 0, (2.14) 

9 t (nM a )+9 / j(nM a M /5 ) = -9 a /? z ' ;flfl/s + Kna a (V 2 n)+nfl a +va / 3[n(a a M / j + 9 / 5M a )] (2.15) 

with kinematic viscosity v = 9t. The terms — d 0l p waals +Knd cl (V 2 n) at the r.h.s. of the 
Navier-Stokes equation can also be written in the form —daP K a where the pressure tensor 
P x p is related to the free energy functional of the van der Waals fluid [30] 

Y = Jdr[ip(n,9) + ^(Vn) 2 ], (2.16) 

ip(n,6) being the bulk free energy density 

3n 9 

xp(n,9)=ne]n(- )--n 2 . (2.17) 

o — n o 

The pressure tensor is [31] 

Pxp = V$<ip+KdKndpn (2.18) 

with 

p = p waals - K nV 2 n-^(Vn) 2 (2.19) 

where p waals = nip (n) — xpis the equation of state with the critical point at n c = 1, 9 C = 1. 

In the sequel, we will consider the case of an acceleration directed upwards: a x = 
0, ciy = g. Periodical boundary conditions were considered in the horizontal direction 
and standard bounce back boundary conditions [8-12] were imposed on top and bottom 
walls. 



3 Numerical results 

In this Section we report the results of our simulations. For runs we used either a square 
lattice with A x = Ay = 1024 or a rectangular one with A x = 512 and Ay = 4096, lattice 
spacing Ss = 1/256 and time step St = 10~ 5 , as in [7]. All quenches below the critical 
temperature 9 C = 1 were to the temperature 9 = 0.79 where the coexisting densities are 
n liquid = 1-956 and n vavor = 0.226. Each simulation was started with small fluctuations 
(0.1%) in the density about a mean value n that was either symmetric (n = 1.09, liquid 
fraction /3 = 0.5) or slightly off-symmetric (n = 1.0, /3 = 0.45). The parameter k controlling 
the surface tension was set to 5 x 10~ 6 to have an interface thickness of ~6 lattice spacings. 
The corresponding value of the surface tension a was evaluated by using the Laplace law 
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[25] and a value cr = 2.0 x 10~ 3 was measured. The value of the constant g controlling the 
external acceleration, was varied in the range [0.0001,0.005]. Results were not dependent 
on its specific value as shown in the following. The viscosity was varied by changing 
t. We used the values r = 10 -4 and T = 10~ 3 , which allow us to access low and high 
viscosity regimes, respectively, as shown in a previous work where the present model 
without gravity was used to study the phase separation in liquid-gas systems [7]. 

The process of phase separation depends on the interplay among three main driving 
forces: The viscous one F v = hvIrCr, the gravity one F g = ngl\, and the surface tension 
one F s = ctIr. It may be then useful to evaluate their relative contributions introducing the 

Bond number Bo = — = — R , the Capillary number Ca = — ^ = — — , and the ratio — = 

F s a F s a Ca 

F gl 2 

-f = — — . The choice of the input parameters is such to access the ranges 0.05 <Bo< 2.5, 
F v vcr 

0.04 < Ca < 0.4, and 0.13 <Bo/Ca< 65. 

After the initial stages when the mixture starts to phase separate, the effect of the 
gravitational force is to accumulate material at walls: The heavy phase (liquid) is moved 
to the top wall and the light phase (vapor) stays at the bottom wall. The evolution of 
domains in the cases at high (t = 10~ 3 ) and low (t=10~ 4 ) viscosities is shown in Figs.QJS] 
for g = 0.005. At intermediate times between t ~ 3 and t ~ 10, anisotropic patterns can be 
observed in the bulk region far from the walls, with domains slightly elongated along the 
vertical direction. The main difference that can be observed between the two figures is 
the presence of many droplets in the case at high viscosity as compared to the case at low 
viscosity. The reason is due to the fact that in the latter case hydrodynamics is effective 
in coalescing droplets, thus producing a more homogeneous pattern. 

In order to gain some insight into the law governing the accumulation of material 
at walls, we measured the average thickness L of layers adjacent to the walls. For each 
column of the lattice we looked for all the sites next to the walls where there was an 
interface between the liquid and vapor phases. To be more specific we found all the lattice 
sites along the x-direction with the smallest distance yt (x) to the bottom wall such that 
[n{x,yl(x))—n][n(x,y* h (x) +1) — n] <0 with y\ (x) < A y 1 2 and all the sites with the smallest 
distance A y — y* t (x) to the top wall such that [n(x,y* t (x)) — n][n(x,y* t (x) — 1) — ft] <0 with 
y*{x) > Ay/2. We defined 

1 A A y 



L = 2A E M W + Cfy - y? (*))]< -f- (3-D 

x x=l 

The time evolution of L is reported in Fig. [3] for different values of r and g. We have 
a clear and convincing indication that in all the cases the growth is consistent with a 
power law with growth exponent a = 1 which depends neither on g nor on x. The growth 
is observed over almost two time decades until the system is entirely separated in two 
parts of different composition. In particular, we want to stress the fact that the exponent 
oc = 1 is observed in both cases, when the gravitational force Fg is small compared to the 
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Figure 3: Time evolution of the average thickness of layers at walls in the cases with g = 
0.0001(A),0.001(D),0.005(o) for r = 1CT 3 (filled symbols), 1CT 4 (empty symbols) and £ = 0.5. The solid 
line is a guide to the eye and has slope 1. 



surface tension force F s (Bo < 1), as well as when F ? is small compared to the viscous 
one F v (Bo/Ca < 1). This indicates that the existence of the same scaling exponent a 
in the gravity direction is exclusively due to the presence of the gravity force F g > 0. 
Our result is in agreement with previous studies of mixtures where hydrodynamics was 
neglected [6,32-34]. In another study on the phase separation of binary fluids [35] where 
hydrodynamic effects were considered, it was argued that the growth exponent is a = 
0.6±0.1 and is not affected by the presence of gravity. The present study shows that 
the value of the growth exponent is independent on the value of the viscosity and of 
the gravity. Similar results were obtained when considering the case of a slightly off- 
symmetric mixtures with ft = 0.45. 

In order to better characterize the morphology of domains we simulated the be- 
haviour of a very large system with symmetric composition (j6 = 0.5) and size 512 x 4096, 
for g = 0.005 and r £ {10~ 4 ,10~ 3 }. Also in this case we found that L grows with the 
exponent 1. In order to estimate the domains size in the two spatial directions, due to 
the anisotropy induced by gravity, we computed the inverse of the first moment of the 
structure factor [36, 37] 



(3.2) 
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Figure 4: Time evolution of the size of domains Ry (upper panel) and R x (lower panel) in the case with g=0.005 
for t=10~ 3 (•), 10~ 4 (o) and /3 = 0.5 on a lattice 512x4096. The lines serve as eyeguide and have the slopes 
1/2 (dotted line), 2/3 (dashed line), and 1 (full line). 
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and similarly for the y direction where 

C(k,f)=<n(k,t)n(-k,t)> (3.3) 

with n(k,t) the spatial Fourier transform of n{x,t) — h. The results are shown in Fig. |4j 
Along the gravity direction, we see that Ry starts to grow with the same exponents as 
in the case without gravity: 1/2 for t = 0.001 and 2/3 for t = 0.0001 [7]. The transition 
to a regime consistent with the growth exponent 1 is observed at later times. Along the 
horizontal direction, we get R x ~ £ 1/<2 in the high viscosity case (t = 0.001), a behaviour 
that is similar to the case when the system is subjected to no gravity. However, at low 
viscosity (t = 0.0001) the growth exponent along the horizontal direction in the presence 
of gravity is smaller than the expected value 2/3, which is achieved without gravity. It 
is interesting to note that similar conclusions were drawn in the case of phase separa- 
tion of binary mixtures under gravity [38]. In their study, the authors observed that in 
the diffusive and viscous regimes the growth exponent is always equal to 1 along the 
vertical (gravity) direction, while in the horizontal direction there is a slowing down of 
the growth rate with respect to the case without gravity [38]. In order to better elucidate 
these features one would need to perform higher resolution simulations to access a wider 
range of length scales. 

4 Conclusions 

In this paper we have introduced an external gravitational force in an isothermal lattice 
Boltzmann model for the van der Waals fluid. We have studied phase separation in sys- 
tems with different viscosities and various values of the gravitational acceleration. In 
the absence of gravity, the growth exponent is known to have specific values [7], which 
depend on the fluid viscosity. When the liquid-vapor system was subjected to the grav- 
itational force, we measured the evolution of the characteristic size (along the gravity 
direction) of the growing domains and found the same exponent a = 1 for all the cases 
considered, even if the fluid viscosity and the gravitational acceleration were different. 
Further extension of our parallel computing code to three dimensions would allow to 
evaluate the growth exponents in a more realistic case. 
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